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ABSTRACT 

We develop and analyze a simple cellular automaton (CA) model that re- 
produces the main properties of the evolution of soft X-ray coronal loops. We 
are motivated by the observation that these loops evolve in three distinguishable 
phases that suggest the development, maintainanee, and decay of a self-organized 
system. The model is based on the idea that loops are made of elemental strands 
that are heated by the relaxation of magnetic stress in the form of nanoflares. 

In this vision, usually called “the Parker conjecture” (Parker 1988), the origin 
of stress is the displacement of the strand footpoints due to photospheric con- 
vective motions. Modeling the response and evolution of the plasma we obtain 
synthetic light curves that have the same characteristic properties (intensity, 
fluctuations, and timescales) as the observed cases. We study the dependence of 
these properties on the model parameters and find scaling laws that can be used 
as observational predictions of the model. We discuss the implications of our 
results for the interpretation of recent loop observations in different wavelengths. 

Subject headings: Sun: corona - Sun: flares - Sun: magnetic topology - Sun: 
X-rays, gamma rays 


1. Introduction 

EUV and soft X-ray observations of the solar corona reveal a high degree of structuring 
and intermittency, especially in and around active regions (ARs). Due to the frozen- in 
condition of the coronal plasma, the emitting material is ordered by the magnetic field 
in elongated features commonly known as coronal loops. The study of the structure and 


1 Institute de Astronomfa y Flsica del Espacio, CONICET-UBA, CC. 67, Sue. 28. 1428 Buenos Aires, 
Argentina 

2 Member of the Carrera del Investigador Cientifico y Tecnologico, Consejo Nacional de Investigaciones 
Cientificas y Tecnicas (CONICET), Argentina 

3 NASA Goddard Space Flight Center, Code 671, Greenbelt, MD 20771, USA 



2 


evolution of loops is therefore fundamental for understanding coronal dynamics, particularly 
in relation to the origin of coronal heating (Mandrini et al. 2000). Historically, the first 
attempts to physically describe coronal loops began with what is known as the Rosner, Tucker 
& Vaianna (RTV, 1978) approach, which considers loops as monolithic structures containing 
plasma in static equilibrium. Although soft X-ray loop observations seem to be consistent 
with the RTV assumptions (Porter & Klimchuk 1995), more recent observations from the 
Transition Region and Coronal Explorer (TRACE) suggest that EUV loops are actually too 
dense to be in static or quasi-static equilibrium (Aschwanden et al. 2001, Winebarger et 
al. 2003). The properties of these loops are more consistent with the predictions of models 
based on impulsive heating (Winebarger et al. 2003, Klimchuk 2006). 

TRACE observations also show a high degree of filamentation of the coronal structure, 
suggesting that what appear to be monolithic loops can actually be collections of unre- 
solved individual strands with more or less independent evolutions. Although the question 
of whether coronal loops are isothermal or multithermal has been a matter of a heated 
debate during several years (Aschwanden & Nightingale 2005, Schmelz & Martens 2006), 
recent observations seem to indicate that both kinds of loops actually occur (Schmelz et al. 
2009). Combining all the available information, there are strong arguments in favor of many 
loops being made of unresolved strands that evolve more dynamically than what quasi-static 
models predict (for a recent review on the subject see Klimchuk 2009). 

There are two main families of mechanisms proposed to explain the origin of coronal 
heating: the so called AC and DC heatings. Generally speaking, AC refers to the dissipation 
of MHD waves, while DC relates to the dissipation of highly stressed magnetic fields. The 
latter involve the release of free magnetic energy at locations of strong gradients in the 
magnetic field, i.e. , current sheets. Ultimately, the source of energy is the distortion of the 
magnetic structures by convective motions. This is precisely the mechanism envisioned by 
Parker (1988). He proposed that, since the footpoints of magnetic strands are being displaced 
by photospheric granular motions, there must be a continuous increase of the magnetic stress 
between neighbor strands. These stresses accumulate until some threshold is reached and 
reconnection occurs. In the process, stored magnetic energy is released and the plasma is 
heated. Because of similarities with the usual flaring process and the energy scale involved, 
these reconnection events are usually called nanoflares (Cargill 1994, Cargill & Klimchuk 
2004) . One key point of the process described above is the existence of a threshold or critical 
value for the reconnection to occur. A theoretical support for the existence of this threshold 
is given in Dahlburg et al. (2005, 2009). They found that when certain field misalignment is 
reached the secondary instability provides the mechanism for a rapid energy release. 


The intermittent and filamentary behavior, and the possible presence of power law 
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distributions (see e.g., Aschwanden k Parnell 2002), led some authors to propose that the 
solar corona might be a good example of a self-organized critical (SOC) system (Bak 1990). 
Beginning with the seminal work of Lu k Hamilton (1991), there have been an important 
number of studies involving the use of SOC models in relation to the coronal heating problem 
(for a review on SOC models see Charbonneau et al. 2001). 

In our recent paper, Lopez Fuentes et al. (2007), we studied the evolution of coronal 
loops observed with the Solar X-ray Imager (SXI) on board the Geosynchronous Opera- 
tional Environmental Satellite 12 (GOES- 12). Since this instrument observes the sun almost 
continuously, it allowed us to follow the evolution of a set of coronal loops from birth to disap- 
pearance. We found that the loop evolution can be separated into three main phases, which 
we called: rise, main and decay. As we indicated there and began to explore in Klimchuk, 
Lopez Fuentes k DeVore (2006), the observed evolution can be related to the development, 
maintenance and destruction of a critical system. Here, we develop a simple model based on 
a cellular automaton (CA) that reproduces the observed evolution. We model the plasma 
response to the heating using the EBTEL hydrodynamics code (Klimchuk et al. 2008). We 
then create synthetic light curves taking into account the temperature-dependent response of 
the SXI instrument. We compare these with the observed light curves and analyze how the 
different parameters of the model influence the properties of the obtained light curves, such 
as the characteristic intensities and timescales, which we studied in detail in Lopez Fuentes 
et al. (2007) for the observed cases. 

In Section 2 we sumarize the results from Lopez Fuentes et al. (2007) and show examples 
of observed light curves that will be later compared with the model. In Section 3 we describe 
the CA model and in Section 4 we analyze our results. We discuss and conclude in Section 5. 


2. Observed loop evolution 

In Lopez Fuentes et al. (2007) we studied the evolution of a set of coronal loops observed 
with the Solar X-ray Imager (SXI) on board GOES-12 (Hill et al. 2005). This instrument 
has the advantage of observing the sun almost continuously, allowing us to follow the whole 
evolution of loops. Previous similar studies based on Yohkoh/SXT or TRACE observations 
were limited by the occultation times of the satellites and times of interrupted observation due 
to high radiation levels (e.g., when passing through the South Atlantic anomaly), therefore 
reducing the continuous observation times to a small portion of the loop lifetime (see e.g. 
Porter k Klimchuk 1995). 

We analyzed a set of 457 images taken with the telescope in the filter position “open” , 
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which is most sensitive to emission from plasmas below 2MK. The analyzed data span three 
days (August 26-28. 2003) of observations with a cadence of approximately 6 images per 
hour. Prom these data we selected and tracked the evolution of 7 coronal loops. Following 
a procedure that is thoroughly described in Lopez Fuentes et al. (2007), we obtained light 
curves of the intrinsic (background subtracted) intensity of each of these loops. We found 
that the loop evolution can be separated into three parts: a rise phase during which the loop 
intensity slowly increases, a main phase during which the intensity remains approximately 
constant (or its long term variation is much less pronounced), and a decay phase during 
which the intensity decreases until the loop is almost indistinguishable from the background. 
In Figure 1 (upper panels) we reproduce adapted versions of two of the light curves studied 
in Lopez Fuentes et al. (2007). In the original analysis we applied a procedure to remove the 
background contribution that left a residual component of long term variation (see discussion 
in Section 3.1 in Lopez Fuentes et al. 2007). For an easier comparison with the model light 
curves, here we removed the residual background by subtracting in each case a linear function 
of time, so that the initial and final intensities of the loops are reduced nearly to zero. 

We found that the durations and characteristic timescales of the rise and decay phases 
are longer than the expected cooling times. We concluded that this is consistent with two 
alternate scenarios. In the first one, loops are monolithic and the heating source varies very 
slowly with time, so the plasma responds quasi-statically. In this way, the three evolution 
phases would reflect the phases of the long term change of a slowly varying heating source. 
In the second scenario, loops are made of unresolved and independent magnetic strands 
which are heated by impulsive events such as nanoflares. Since the observed loop emission is 
the summed contribution of the constituent strands, the different phases would be directly 
related to the variation of the number and intensity of the energy dissipation events that 
heat the strands. We argued that under this particular scenario the observed loop evolution 
can be related to the development, maintenance and destruction of a self-organized system. 
In the next sections we present and analyze a model that reproduces precisely this kind of 
behavior. 


3. Cellular automaton model 
3.1. General description 

In this Section we describe the construction of a cellular automaton (CA) model that 
reproduces the main features of the observed loop evolution. Motivated by the arguments 
given in Section 1, the model simulates the process described there regarding the scheme: 
motion of magnetic strand footpoints — + increase of magnetic stress — > energy release by 
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reconnection. Figure 2 shows a cartoon that illustrates the main idea. 

Let us begin with a bunch of thin parallel magnetic flux tubes or strands that connect 
two different regions of the solar surface (panel (a)). For simplicity the strands are drawn 
straight, so the top and bottom planes correspond to two relatively distant photospheric 
regions of opposite magnetic polarity. As time goes on, strand footpoints are horizontally 
dragged by random photospheric motions, so the strands become progressively more mis- 
aligned. This produces growing magnetic stresses and intensifying current sheets in the region 
separating adjacent strands (panel (b)). For two given neighbor strands, the misalignment 
angle increases until a critical value is reached (i.e., the secondary instability develops) and 
the energy is released in the form of an impulsive event. Consequently, the field between the 
two strands relaxes. This sudden change in the field may occasionally cause new instability 
conditions between these strands and their other neighbors. Under some conditions instabil- 
ities propagate farther producing avalanches of events involving many strands. This is the 
typical situation considered in SOC models. For clarity, we adopt the following terminology. 
The basic interaction between two strands is called a “reconnection event.” The total release 
of energy into a single strand by its reconnection with one or more neighbors over a short 
period of time is called a “nanoflare.” This follows the convention used by the loop modeling 
community. 

The system continues to evolve in response to the footpoint driving until there is a 
rather abrupt appearance of a critical state in which the energy input by the photospheric 
motions is statistically balanced by the dissipation due to reconnection events. In the critical 
state the total energy of the system fluctuates intermittently around a constant level. Local 
intermittency results from the occurrence of nanoflares of different sizes. 

The idea is then to develop a numerical algorithm that quantitatively reproduces the 
above scheme. First of all, we need to translate footpoint displacements to magnetic field 
stresses. To accomplish this we begin with the simple scenario shown in Figure 2 (c). Let us 
concentrate in one particular strand and, for simplicity, suppose that only one footpoint (the 
lower one) is being displaced. If the particular strand has a vertical magnetic field strength 
B v , any horizontal displacement of its footpoint will produce a new horizontal component of 
the field 


5B h = B v tan 9, (1) 

where 6 is the angle that the now inclined strand forms with the vertical direction. It is easy 
to see that if the strand has a length L and the footpoint has been displaced a distance d, 
then tan 9 = d/L. Therefore: 
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5B h = B v d/L. (2) 

The distance d is a characteristic length scale of the photospheric motions, which can be 
associated with the sustained displacement of a strand footpoint during the turnover time 
of a convective granule (~ 10 3 s). Assuming a typical photospheric velocity v ~ 1 km/sec, 
results ind~ 1000 km. After a number of successive displacements the strand would acquire 
a horizontal field 

B h = T,SB h . (3) 

At this point we have assumed that Bh increases with each and every step. This is only true 
if the new step is in the same direction as the previous step or if the change in direction 
causes the strand to wrap around a neighbor. The latter will generally be true if the random 
walk step size is larger than the characteristic separation of neighbor footpoints. However, 
even then, some steps will cause the strand to back track and Bh will decrease. We treat 
this more realistic situation below. 

As a consequence of the random walk, the horizontal field Bh will gradually increase 
until the critical condition referred to above is reached. If we call 9 C the critical angle of 
misalignment between two adjacent strands, it is easy to see that there is an associated 
critical value of the difference in their horizontal field components: 

B c = B v tanO c . (4) 

Note that 9 C refers to the relative angle between the adjacent strands and not the inclination 
from vertical, 9 , in Equation 1. 

Assuming ID displacements (see Section 3.2), we relate the misalignment between two 
neighbor strands (with indices p and q) to the difference of their horizontal field components 
A B h = B h (p ) — B h (q). The fulfillment of the critical condition implies | A5^| > B c . If we 
assume that all the field is redistributed after reconnection, the new (primed) horizontal field 
in each strand will respectively be: 

B' h {p) = B h (p ) - AB h /2 and B' h {q) = B h (q) + AB h /2. (5) 

The released energy corresponds then to the difference in the total magnetic energy 
before and after reconnection. It can be easily demonstrated that 



7 


A E = (. E' 


E) 


1 

8tt 


(A B h f 


( 6 ) 


per unit volume, where we have assumed that the energy is divided equally between the 
plasma in the two strands. 


3.2. Model implementation 

For the numerical implementation of the model we use a ID array of N m points. Each 
point corresponds to a single magnetic strand and, initially, all strands have horizontal 
magnetic field Bh — 0. At each time step, a small random amount of field 5B h is added to 
each site according to a statistical distribution that simulates the photospheric displacements 
of strand footpoints. Although the field increase in the model is ID, we apply a distribution 
based on 2D footpoint displacements. The idea is illustrated in Figure 3-a. If one compares 
the position of a strand footpoint at time step (t — 1) with its position at time step t, the 
probability that during the following time step, (t + 1) the footpoint moves farther away 
increasing the horizontal magnetic strength is approximately 3/4, while the probability that 
it moves back on its previous track, canceling (or almost canceling) the previous displacement 
is roughly 1/4. Therefore, we add the 5B h values according to the distribution shown in 
Figure 3-b. Values from the gaussian distribution centered on d 0 are 3 times more probable 
than values given by the gaussian centered on —d 0 . The distance do = 1000 km is the average 
horizontal displacement of a strand footpoint during the turnover time of a granular cell (in 
our scheme, the duration of each time step). The width of the gaussians is 20% of d 0 . As 
explained in the previous Subsection, the random displacements taken from this distribution 
translate into magnetic inputs through Equation 2. 

It is easy to see that under this simple ID scheme all strands will eventually acquire 
positive values of B h . To keep the simplicity of the model and to be consistent with the fact 
that neighbor strand footpoints tend to move in different directions, we alternate the sign of 
the displacements d among neighbor strands. This means that, identifying each mesh site 
with an index i, we add the random values obtained from the distribution to strands with 
even i, and we subtract them from strands with odd i. Thus, in the long run even strands 
acquire positive B h while odd strands acquire negative B h . This alternation of 5B h simulates 
in ID the magnetic stress produced by footpoints moving in different random directions. 

The efficiency of the process that injects “horizontal component” of the magnetic field 
through the motion of the footpoints depends on the average distance between neighboring 
strands. If the strands are closer than d 0 , then the process can be considered efficient in 
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tangling the strands and accumulating magnetic stress. However, if there is no injection of 
new flux, then the distance between strands will increase as the footpoints disperse. The 
driving process will become progressively less efficient. To simulate the decay of the driving 
mechanism due to footpoint dispersion, at a chosen point during the system evolution we 
make the distribution shown in Figure 3-b slowly evolve to the symmetrical distribution 
shown in Figure 3-c. In the final state, the driver distribution is completely inefficient at 
creating new B h . As we will see below, this evolution of the driver leads to the decay of 
the system. The start time and the duration of the evolution from the initial to the final 
distribution are set as parameters of the model. 

In Subsection 3.1 we defined B c = B v tan 6b . Therefore, tan# c is one of the relevant 
parameters of the model. At each time step, after adding the random SBh values to all the 
strands, we compare Bh on each strand with its immediate neighbors. For a strand having 
index i we compute the difference A B h (i) = B h {i) — B h (i + 1). We do this for all the strands 
in the array. We use periodic boundary conditions, so strand i = N m is compared with 
strand i = 1, and in that case ABh(N m ) = B h (N m ) — B h (l). Next, we identify all the indices 
where \AB h (i)\ > B c . For all the sites where this happens, we redistribute the horizontal 
field according to Equation 5. 

As we discussed in Subsection 3.1, each field redistribution is associated with a recon- 
nection event that frees an energy A E, in erg cm" 3 , according to Equation 6. This energy 
is used to heat the plasma in both strands involved in the event. Before proceeding to the 
next time step we repeat the previous test and corresponding field redistribution as many 
times as necessary until all j Ai?^(z) | values are less than B c . Once this state is reached, the 
system is allowed to evolve to the next time step, a new set of random SBh values is added 
to the Bh of the strands, and the process described above is repeated. 

During the evolution of the system we keep record of the energy dissipated in each 
strand at each time step. Here, we assume that the reconnection events (A E) occur on a 
timescale that is much shorter than the time step duration ( 10 3 s). As we defined in 

Section 3.1, the sum of all the A E “energy packages” that heat the plasma in the strand 
during a given time step is called a nanoflare. The final output of the CA model is the array 
of nanoflare energies dissipated in each strand as a function of time. 

In the above description we assumed that when the critical condition (\AB h \ > B c ) 
is fulfilled all the available field AB/, is distributed (reconnected) among the two strands. 
This is not realistic. On the Sun, each strand of the tangled field is wrapped around many 
other strands. It is only the horizontal field component along the section of mutual contact 
that is free to reconnect and heat the plasma. We account for this in our simple model 
by introducing a parameter /. The initial response to the reconnection is to change the 
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horizontal field locally by an amount AB h /2. This change is subsequently spread out along 
the entire length of strand. The adjustment happens at the Alfven speed and is very rapid. 
The end result is that the horizontal field changes by fAB h /2 along the entire strand. It 
can be demonstrated that in this case the energy release on each of the two strands per unit 
volume is 


A£ =(sM 

Note that the parameter / affects the time required for a pair of adjacent strands to 
once again reach the critical condition through continued driving. If / is small, the strands 
will never be far from critical and nanoflares will recur frequently. If / is large, it will take 
many steps to once again reach the critical level of stress. Of course the nanoflares will be 
weaker in the former case than in the latter, so the temporally averaged heating will not 
be greatly different. There are small differences though, reflected in the / dependence of 
equation (11) discussed later. 

For the analysis, we vary the values of the model parameters B v , L, tan# c , /, and N m . 
Table 1 shows the different values used. For the vertical magnetic field we use values in 
the range 20-120 Gauss, typical of AR loops (Mandrini et al. 2000). We vary the length of 
the loops also considering typical AR values (20-120 Mm). The tangent of the critical angle 
(tan# c ) varies between 0.1 and 1.0, corresponding approximately to angles: 6 deg, 11 deg, 
22 deg, 31 deg, 39 deg, and 45 deg. The reconnection factor / varies between 0.1 and 1.0, 
meaning that 10% to 100% of the available magnetic stress is released by reconnection. The 
total number of strands, N m , goes from 20 to 1000, which is roughly the number of elemental 
magnetic flux tubes believed to be present in a coronal loop. We vary one parameter at a 
time keeping the rest of them at the intermediate values shown in bold face in Table 1. 

As we stated above, the output of the CA model is an array of nanoflare energies dissi- 
pated in each strand as a function of time, with time steps of ~ 10 3 s. For any combination 
of model parameters the initial evolution of each strand is such that the number of reconnec- 
tion events increases from zero (rise phase) to a stationary level (main phase). The system 
reaches this state abruptly, not asymptotically, and maintains it for as long as the driver 
follows the distribution of Figure 3-b. As soon as the driver begins to evolve towards the 
distribution of panel c, the number of events decays until they almost disappear at the end 
of the system’s evolution. Therefore, the evolution of the full system qualitatively follows 
the evolution of the observed coronal loops in terms of rise, main, and decay phases. 


A (AB0 : 


( 7 ) 



3.3. Plasma response and light curve construction 


Since our main motivation is to compare the output of the CA model with the loops 
studied in Lopez Puentes et al. (2007). the next step is to transform this output into an 
observed signal such as the light curve of a coronal loop. To do this, we simulate the plasma 
response using the EBTEL (Enthalpy-Based Thermal Evolution of Loops) model, developed 
by Klimchuk, Patsourakos & Cargill (2008). The temperature and density of the coronal 
plasma tend to be reasonably uniform along the magnetic field, and EBTEL computes the 
evolution of the spatially-averaged values. This approximate 0D solution is quite similar to 
the exact solution given by ID hydrodynamic codes, but it takes 3-4 orders of magnitude 
less time to compute. EBTEL is therefore ideally suited for our study in which we simulate 
the evolution of thousands of individual strands. 

The main input of EBTEL is the heating rate profile, i.e., the variation of the heating 
rate with time in the strand. Here, we assume that nanoflares have a triangular profile of 
a given duration r. The upper panel of Figure 4 shows a series of “triangular” nanoflares 
that heat one of the strands of the CA model during a portion of the main phase. The time 
integral of each triangular heat function corresponds to the total energy provided by the 
particular nanoflare. The duration of the nanoflares shown in Figure 4 is 1000 s. This means 
that, in this case, the available heat is distributed (with a triangular profile) along a full 
time step of the CA model. The other two panels of Figure 4 show the output of EBTEL: 
the plasma temperature and density. During the time interval covered in the panels there 
are both isolated events (like the nanoflare starting at t = 4.3 x 10 4 s), and “trains” of 
events occurring in rapid succession (such as the one starting at t = 4.7 x 10 4 s). Notice 
that while the plasma relaxes almost fully after the isolated event, during the succession 
of multiple events the plasma is still in a state of high density and temperature due to the 
previous nanoflare when the next event begins. As we will see in the following sections, this 
has consequences in the way the plasma emits and the observed loop intensity evolves. 

The loops studied in Lopez Fuentes et al. (2007) have measured densities in the range 
0.9-1.9xl0 9 cm -3 , and measured temperatures between 1.2 and 2.1 MK. These values are 
reproduced by parameter values in the ranges given in Table 1. 

Once we have the temperature and density evolution of each strand, we compute the 
observed intensity using the emission measure and the known temperature dependent sen- 
sitivity of the instrument. The emission measure per pixel of a single unresolved strand is 
EM = n 2 AV, where AV is the volume of the strand portion covered by a pixel. This is 
roughly AV = l P i X A s , where l P i X is the pixel dimension and A s is the cross sectional area 
of the strand. The contribution to the intensity in a pixel is then I = EM S(T), where 
S(T) is the SXI response function. We sum the contributions of all the strands to determine 
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the intensity of the loop and. finally, we multiply this intensity by a factor that accounts 
for the proportion of the loop covered by a single pixel. The final output is the observed 
intensity of the loop per pixel per second as a function of time. The data are sampled using 
the instrument temporal resolution to simulate observed light curves. 


4. Results 

The lower panels of Figure 1 show two model light curves obtained in the way described 
in Section 3.3. The cases shown have the same properties as the observed light curves of 
the upper panels. Given the random nature of the nanoflares, we do not expect an exact fit 
of the model to the observed data. For the comparison we instead consider the light-curve 
properties studied in Lopez Fuentes et al. (2007) (see Table 2). 

The values of the model parameters used to construct the synthetic light curves of 
Figure 1 are, for Loop 1: B v = 37 G, tan# c = 0.15 and / = 0.2; and for Loop 4: B v = 20 G, 
tan# c = 0.2 and / = 0.1. For both cases we use N m = 100, and r = 1000 s. For the loop 
length and diameter we use the measures obtained in Lopez Fuentes et al. (2007): L = 110 
Mm for loop 1, L = 70 Mm for loop 4, and a diameter of 1.8 x 10 9 cm for both. 

We compute for both, observed and synthetic light curves, the mean intensity of the 
main phase ( I m ) and the duration and timescales of the rise and decay phases. To obtain 
these values we first identify, by visual inspection, the approximate start and end times of 
each phase. Using the same definitions as in Lopez Fuentes et al. (2007), we have for the 
rise phase timescale: 


rise- (A//AW)' w 

Here, I is the mean intensity during the rise phase, A I is the intensity variation and A t r i se is 
the duration. A similar definition is used for the decay phase. As a proxy for the amplitude 
of the intensity fluctuations, we compute the root-mean-square (RMS) of the intensity 
variation during the main phase relative to its ten-point running average. In Table 2 we 
compare these properties for observed and modeled light curves. It can be seen that the 
model reproduces the light curve properties consistently. 

We remind the reader that the observed light curves shown in Figure 1 have been further 
processed (residual background removed, see Section 2) from the versions studied in Lopez 
Fuentes et al. (2007). Therefore, the properties listed in Table 2 do not necessarily coincide 
with the corresponding values presented in the previous paper. 
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Regarding the durations and timescales it is interesting to see that, since the light 
curves have initial intensities very close to zero. I in the Equation 8 is approximately A//2. 
Therefore, 

rp ^rise / q\ 

-‘-rise ~ 2 ’ V'-v 

We then expect the timescales to be approximately one half of the durations (compare the 
corresponding values in Table 2 for both the rise and decay phases). The relation from 
Equation 9 is not valid, in general, for the light curves studied in Lopez Fuentes et al. 
(2007), in which there is still a remnant background component of the intensity. 

To investigate how the model parameters determine the properties of the loops, we vary 
each of them separately according to Table 1 while keeping the rest fixed at the values shown 
in bold face (see Section 3). In Figure 5 we show log- log plots of the the mean energy per 
unit volume released in reconnection events, A E (see Equation 7), versus different model 
parameters. The lines are least-square fits of the data, and their slopes are indicated in the 
panels. According to these plots, in our CA model the energy of the reconnection events 
approximately scales with B% m and (tan0 c ) 1-89 . This can be easily explained by noting that 
in Equation 7, when reconnection occurs, A B h « B c = B v tan9 c . Therefore, it is expected 
that A E oc Bl (tan# c ) 2 . 

The horizontal field decreases by an amount fAB h /2 during a reconnection event. Sev- 
eral timesteps are typically required to rebuild the field to the critical level whereupon 
another event occurs. From equation (2) we see that the number of recovery timesteps is 
given by 


N r 


— — tan 9 C . 
2 d 


( 10 ) 


From this we obtain the mean nanoflare heating rate (averaged both in time and along the 
strand): 


A E _( f\ vBlt&n9 c 

N r t s ~ V V L 


( 11 ) 


where t s is the time step duration and v = d/t s is the photospheric driving velocity. We have 
assumed that there is one reconnection event per nanoflare, whereas often there are two, but 
this is not important for identifying scaling relationships. In Figure 6 we show log- log plots 
of the mean heating rate versus B v , tan 9 C , and L. The dependence on the parameters is: 
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B 202 , L L1 , and (tan# c )°- 97 , which is clearly explained by Equation 11. It is easy to see that 
(Q) does not depend on the number of strands, N m , since both the driver and the critical 
condition are applied to the strands individually. 

We can also analyze scaling relations for the properties of the light curves obtained 
with the model. For each combination of the parameters from Table 1 we obtain synthetic 
light curves following the procedure described in Section 3.3. In Figure 7 we present log-log 
plots of the mean intensity of the main phase, I m , as a function of the same parameters 
as in Figure 6. The scalings found in this case are: B 2 - 37 , L -a8 , and (tan0 c ) L15 . Clearly, 
I m is related to the nanoflare energy, but the dependence is complex and is affected by 
both the detailed response of the plasma and the nonuniform temperature sensitivity of the 
instrument. I m is not linearly proportional to the total radiation emitted by the strand. 
To understand the origin of the dependence, let us now consider the case of quasi-static 
equilibrium. This is a reasonable representation if the time interval between nanoflares is 
short compared to a cooling time. When such conditions apply, the conductive and the 
radiation terms of the energy balance equation are comparable in magnitude to each other 
and to the mean heating rate, (Q) (Vesecky et al. 1979): 

2 T 7/2 

{ Q ) ~ ~ n 2 A 0 T b . (12) 

Here, k 0 is the coefficient of thermal conductivity, and A 0 and b are constants. From the 
first part of Equation 12 we can obtain for the temperature: 

T oc ( Q) 2/ 7 L 4 ' 7 oc B 4 v /7 L 2 / 7 ( tan Q c ) 2 ' 7 . (13) 

Using that b « —0.5 for the temperature range of the SXI loops studied in Lopez Fuentes 
et al. (2007), from equations 12 and 13 we have: 


n oc — oc B 3 J‘ L 3 / 7 (tan 0 C ) 4//7 . (14) 

The observed loop intensity is I = n 2 S(T), where S(T ) is the response function of 
the instrument, which we can express as S(T) = S 0 T a . We found that a « 0.7 for the 
temperature range of the loops studied in Lopez Fuentes et al. (2007). Replacing equations 13 
and 14 in the expression for the loop intensity, we obtain finally: 


Im oc R 2 - 69 L-°- 66 (tan0 c ) 1 ' 34 . 


(15) 
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This expression is consistent with the scalings in Figure 7. Differences are to be expected 
since the strands are sometimes far from quasi-static equilibrium (Fig. 4). 

It might be expected that the intensity of the main phase depends on N m , but this is not 
the case. The reason is as follows. The intensity contribution of each strand is proportional 
to the volume of the strand subtended by the instrument pixel, which is proportional to the 
strand cross-section. Since we compute the strand cross section as the observed cross-section 
of the loop (1.8 x 10 9 cm) divided by N m . an increase of N m means a proportional decrease 
of each strand intensity, so the summed intensity of all strands remains unchanged. 

As we mentioned above, another important property of the light curves is the amplitude 
of the intensity fluctuations, which we quantify from the RMS of the intensity variation 
during the main phase. These fluctuations are related to the frequency and relative sizes of 
the nanoflares. As we briefly discussed in Section 3 the factor /, which relates to the fraction 
of the strand length that is involved in the reconnection, determines the waiting time for 
two neighbor strands to recover the critical condition. If / is small the waiting time will be 
short and the nanoflares will tend to be smaller and more frequent. On the other hand, as / 
tends to 1 the nanoflares will be larger and less frequent. Therefore, it is expected that the 
intensity fluctuations increase with /. The upper panel of Figure 8 shows precisely that. 

The intensity fluctuations are also affected by the number of strands, N m . Since the 
intensity of the loop is the summed contribution of all the strands, given the random nature 
of the release events, the larger the number of strands, the smoother the signal. This is 
shown in Figure 8, middle panel. We see that the RMS intensity fluctuation decreases as 
the number of strands increases. The effect levels off once a sufficiently large number of 
strands are present. 

The RMS fluctuation also decreases as the duration of the nanoflares, r, increases, as 
shown in the lower panel of Figure 8. Individual strands evolve more slowly as the nanoflares 
get longer, and for long enough nanoflares the strands behave in a quasi-static manner. As 
is the case with N m , the effect levels off once the nanoffares are sufficiently long. 

Finally, although the dependence is much weaker than for the previous parameters, the 
RMS variation tends to increase with all the parameter changes that increase the intensity, 
i.e., the increase of B v and tan# c , and the decrease of L. This can be understood in terms of 
the time it takes the plasma to cool fully after a nanoflare. The cooling time varies inversely 
with the nanoflare energy and directly with the loop length (Cargill & Klimchuk 2004). 
Shorter cooling times result in more bursty emission and greater fluctuation amplitudes. 
The increase of the fluctuation amplitude with the increasing intensity is consistent with 
the observational results obtained recently by Sakamoto et al. (2008, and references therein) 
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and discussed also in Vekstein (2009). They also found that the amplitude of the intensity 
fluctuations of Yohkoh/SXT loops is approximately 8% of the value of the loop intensity. 
We find similar values in our own observations and the results obtained with the CA model 
(see the main phase RM S variation in Table 2) . 

Another important property of the light curves is the rise duration. A t rise , which is 
the time it takes to reach the main phase after the intensity first starts to rise. Note that 
the intensity does not increase immediately when footpoint driving begins (Fig. 1). It takes 
time for the first strands to reach the critical condition, even though they have larger-than- 
average steps. We define the waiting time T w to be the elapsed time between the start of 
the evolution (t = 0) and the main phase beginning. It can be easily seen that T w is simply 
the number of steps required to reach the critical condition, \Bh\ = BJ 2, where the factor 
of 2 is because odd and even strands are tilted in opposite directions. The average \6Bh\ per 
timestep that is injected by the driving is only half the value given by equation (2), because 
one in four steps cause the stresses to decrease (Fig. 3). Noting that B c = B v tan# c , we find 
that the waiting time is given by 


T 

-*■ W 


BJ 2 ^ Lta,nd c 
SB h /2 ~ d 


(16) 


in units of the timestep duration (1000 s). In Figure 9 we show log- log plots of T w versus 
the parameters. The scalings are: L 088 and (tan# c ) 0 ' 85 . This is generally consistent with 
equation (16), and we suggest that small deviations are due to the subjective nature of 
identifying the start of the main phase. Notice that T w is expected to be independent of 
B v , which we have confirmed to be true from our simulations. We emphasize again that the 
waiting time is different from the rise duration. The rise duration is more closely related to 
the spread of the driver distribution (e.g., the Gaussian width in Fig. 3b) than it is to the 
time required for the average strand to reach critical. 

The duration of the decay phase depends on both the form of the final driver distribution 
(Fig. 3b) and the timescale for transitioning from the initial to final distributions. We leave 
the details for a later study. 


5. Discussion 

We develop a pseudo-ID model based on a cellular automaton (CA) that reproduces 
the main properties of the evolution of coronal loops observed in soft X-rays. As we describe 
in Section 3, the model is based on the idea that coronal heating is due to the sudden recon- 
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nection of tangled magnetic strands (Parker 1988). Using typical solar ranges for the model 
parameters we qualitatively reproduce the observed intensity amplitude, fluctuations, and 
evolution timescales of soft X-ray loops studied in Lopez Fuentes et al. (2007). The syn- 
thetic loops obtained with the CA+EBTEL model combination also reproduce the measured 
physical properties of the observed loop plasma such as density and temperature. 

Although our model is related to usual models of self-organized criticality (SOC, see 
e.g. Chabonneau et al. 2000), it differs in some important aspects. In most SOC models the 
rules for the critical condition and the redistribution of the field among neighboring strands 
are based on mean values of the field strength. The strength of each strand is compared 
with the mean value of its neighbors, and if the critical condition is fulfilled the field is 
redistributed equally among those neighbors. One may classify this type of model isotropic. 
As we describe in Section 3, our CA model is anisotropic in this respect. The other important 
difference is that SOC models use null boundary conditions. The mesh points (strands) at 
the boundary are set to have zero magnetic field (what we call Bh in Section 3) at all times. 
Here, we use instead periodic boundary conditions. It may be partly due to these differences 
that the nanoflares produced with our model do not follow power law distributions, which 
are characteristic of SOC models (for an anisotropic SOC model associated with power law 
distributions, see e.g. Morales & Charbonneau 2008). We must stress that our main goal 
here is not to obtain such distributions, but instead to explain the evolution of loops in the 
frame of the nanoflare model. It is worth pointing out, though, that our CA model reaches 
a steady critical state in much the same way as SOC models do. 

As we mention in Section 1, in recent years there has been a debate over the issue 
of whether coronal loops are isothermal or multithermal. The argument goes that if loops 
are isothermal they are likely to be monolithic structures, while multithermality implies that 
loops are composed of different unresolved strands evolving independently. There seems to be 
evidence in favor of both kinds of observations (see e.g., Schmelz et al. 2009), so a consensus 
is forming that both scenarios actually occur. Let us analyze this possibility in terms of the 
model presented here. For instance, if the reconnection mechanism is not very efficient (e.g., 
a small / factor, see Section 3) and all nanoflares are approximately the same size, then 
the heating injection will be difficult to distinguish from a single source of constant heating. 
In that case, the strand will keep more or less the same temperature along the main phase 
of the evolution. If all the strands that comprise the loop have similar physical conditions 
(it is expected that they all have approximately the same field strength, for example), then 
the loop will be nearly isothermal. These loops can be easily taken as examples of quasi- 
static evolution. On the other hand, if nanoflares occur less frequently, so the waiting time 
between consecutive nanoflares is not negligible, the strands have time to cool to lower 
temperatures between one nanoflare and the next. Since nanoflares occur randomly, at any 
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given time during the evolution of the loop there will be strands simultaneously having 
different temperatures. That would correspond to a multithermal loop observation. 

There is yet a third possibility. If the nanoflares occur only once in each strand, with- 
out repeating, and if they all occur at approximately the same time, then the loop will 
appear roughly isothermal at any given moment. The temperature will decrease on a cool- 
ing timescale, however, and the loop will be short lived. We call this scenario a short duration 
nanoflare “storm” (Klimchuk 2009). Many warm EUV loops can be explained in this way. 
Hot soft X-ray loops tend to live much longer than a cooling time, on the other hand (Lopez 
Fuentes et al. 2007). If they are indeed isothermal, then they can only be explained by 
rapidly repeating nanoflares or by truly steady heating. 

The three alternative scenarios described above relate to the question weather the same 
loop is observable by instruments that are sensitive to different temperatures. In the case of 
nanoflares that repeat rapidly in each strand of the loop bundle, the loop will be visible only 
over a narrow range of temperatures. It will be seen either in hot soft X-ray emission or in 
warm EUV emission, but not both. In the case of nanoflares that repeat with a time delay 
longer than the cooling time, the loop will be simultaneously visible in both emissions. In 
the final case of a short duration nanoflare storm, the loop will be visible first in soft X-rays 
and then in the EUV. 

Observational evidence has been presented for all three possibilities. The question of 
which is most prevalent is not yet answered. Loops produced by short duration storms 
have been clearly identified and can be found at many locations within active regions (e.g., 
Ugarte-Urra et al. 2006, 2009). On the other hand, soft X-ray loops are most pronounced in 
the central cores of active regions, where EUV loops are less common (Schmieder et al. 2004). 
There are also reports of soft X-ray and EUV loops that are almost but not exactly co-spatial 
(Nitta 2000, Nagata et al. 2003). We attempted to study the warm emission associated with 
the SXI loops studied in Lopez Fuentes et al. (2007). Unfortunately, TRACE data are not 
available, so we had to rely on lower cadence observations from SOHO/EIT. Since only 
four images per day were taken in each channel, we could only compare the hot and warm 
emission at a few times during the loop evolution. We found that some of the SXI loops 
had an EIT counterpart but others did not. Clearly, more systematic observations of loops 
in different wavelengths are needed. 

It is possible that there are two well differentiated populations of AR loops. Hotter 
and longer-lived loops are preferentially located at the core of ARs (Antiochos et al. 2003, 
Ugarte-Urra et al. 2009), while the periphery is dominated by longer, shorter- lived, cooler 
loops (Del Zanna & Mason 2003). If this is indeed the case in general, then the different 
locations and evolutions of these loop classes may be related to different physical properties 
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at their origin. The model we present in this paper applies specifically to hot loops. Whether 
it can explain warm loops at the periphery of ARs is unclear. We can nonetheless speculate 
on the difference in these two regions. The strength of the coronal magnetic field is generally 
greater in the core of an AR and this will tend to produce stronger heating. Our model 
also depends fundamentally on the properties of footpoint shuffling. It may be that these 
properties vary systematically across the active region. This is an important question that 
we hope will be answered with the latest magnetogram observations from SOT/Hinode. 
These observations will hopefully also confirm our proposal that the slow decay phase of hot 
loops is associated with the dispersal of magnetic footpoints to the point where magnetic 
field tangling is no longer efficient. 

In Section 4 we found a series of scaling laws relating observed properties of the loop 
evolution to different physical parameters of the model. Some of these parameters, such 
as the loop length, can actually be obtained directly from the observations. The strength 
and tilt of the magnetic field at the coronal base can be inferred by combining photospheric 
magnetogram observations with models of the coronal structure (see e.g., Lopez Fuentes et 
al. 2006, Mandrini et al. 2000). Therefore, these scaling laws are predictions of the model 
that can be tested observationally. 

From the discussion in the previous paragraphs we conclude that among future impor- 
tant investigations is the systematic study of loop observations in different wavelengths and 
the analysis of possible scalings of the evolution parameters (timescales, mean intensities, 
fluctuations) with the loop physical properties such as length and magnetic field. From the 
theoretical side, the model presented here is still very simple, and more sophisticated ver- 
sions should be developed, such as 2D models that better account for the vector nature of 
the magnetic field and the fact that each magnetic strand is wrapped around multiple other 
strands. We have already started work in this direction. 


This work was supported by the NASA Supporting Research and Technology and Living 
With a Star Programs. We thank Paul Charbonneau for useful discussions. 
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Table 1: Numerical values used for the parameters of the model. 


B v (Gauss) 

20 

40 

60 

80 

100 

120 

L(Mm ) 

20 

40 

60 

80 

100 

120 

tan# c 

0.1 

0.2 

0.4 

0.6 

0.8 

1.0 

/ 

0.1 

0.2 

0.4 

0.6 

0.8 

1.0 

N m 

20 

50 

100 

200 

500 

800 

T 

100 

500 

1000 

2000 

3500 

5000 


Table 2: Compared properties of the observed and modeled lightcurves shown if Figure 1. 



Loop 1 

Loop 4 


Obs. 

Model 

Obs. 

Model 

Rise phase duration (hr) 

3.35 

3.65 

3.74 

3.16 

Rise phase timescale (hr) 

1.97 

1.88 

2.17 

1.75 

Main phase intensity (DN/pix.sec) 

6.82 

7.59 

2.99 

2.82 

Main phase RMS 

0.10 

0.10 

0.13 

0.09 

Decay phase duration (hr) 

2.98 

3.99 

4.51 

4.19 

Decay phase timescale (hr) 

-1.51 

-2.20 

-2.47 

-2.48 



Intensity (DN/pix.sec) Intensity (DN/plx.sec) 


-22 


Observed loop 1 



Synthetic light curve (loop 1 type) 



Observed loop 4 



Synthetic light curve (loop 4 type) 



Fig. 1. — Evolution of observed and synthetic X-ray loops. The light curves obtained with 
the CA model (lower panels) reproduce the main properties of the evolution of the observed 
loops (upper panels). See Section 4. 
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Fig. 3. — The illustration in the top panel helps to explain the use of the driver distributions 
in the middle and bottom panels. Strand footpoint displacements are obtained from the 
distributions and translated into injection of horizontal magnetic field (see Section 3.2). 
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eters of the model (see Section 4). 
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Fig. 8. — Scatter plots of the root mean square 
versus relevant parameters of the model. 








Fig. 9. — Log-log plots of the waiting 
of the model. 
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